function givens_qr()
% GIVENS_QR Givens旋转QR分解演示
% 
% 演示基于Givens旋转的QR分解算法

fprintf('=== Givens旋转QR分解 ===\n\n');

%% 1. Givens旋转基础
fprintf('1. Givens旋转基础\n');
fprintf('Givens旋转矩阵可以将向量的一个分量置零\n\n');

% 演示基本Givens旋转
x = [3; 4];
fprintf('原始向量 x = [%g; %g]\n', x);
fprintf('向量长度: %.6f\n', norm(x));

% 计算Givens旋转参数
[c, s] = givens_rotation(x(1), x(2));
G = [c, s; -s, c];

fprintf('Givens旋转矩阵 G =\n');
fprintf('[%8.6f  %8.6f]\n', G(1, 1), G(1, 2));
fprintf('[%8.6f  %8.6f]\n', G(2, 1), G(2, 2));

% 应用旋转
y = G * x;
fprintf('旋转后向量 y = [%.6f; %.6f]\n', y);
fprintf('第二个分量应该为零: |y(2)| = %e\n', abs(y(2)));

% 验证旋转矩阵性质
fprintf('G的性质验证:\n');
fprintf('正交性 ||G*G^T - I||: %e\n', norm(G * G' - eye(2), 'fro'));
fprintf('行列式: %.6f (应该为1)\n', det(G));

%% 2. 三维向量的Givens旋转
fprintf('\n2. 三维向量的Givens旋转\n');
fprintf('逐步将向量的分量置零\n\n');

x3 = [1; 2; 3];
fprintf('原始向量 x = [%g; %g; %g]\n', x3);

% 第一步：消除x(3)
[c1, s1] = givens_rotation(x3(2), x3(3));
G1 = [1, 0, 0; 0, c1, s1; 0, -s1, c1];
y1 = G1 * x3;
fprintf('第一步旋转后: y = [%.6f; %.6f; %.6f]\n', y1);

% 第二步：消除y(2)
[c2, s2] = givens_rotation(y1(1), y1(2));
G2 = [c2, s2, 0; -s2, c2, 0; 0, 0, 1];
y2 = G2 * y1;
fprintf('第二步旋转后: y = [%.6f; %.6f; %.6f]\n', y2);
fprintf('最终结果应该是 [±||x||; 0; 0]: [±%.6f; 0; 0]\n', norm(x3));

%% 3. Givens QR分解算法
fprintf('\n3. Givens QR分解算法\n');
fprintf('逐个消除矩阵下三角部分的元素\n\n');

% 测试矩阵
A = [1, 2, 3;
     4, 5, 6;
     7, 8, 9;
     10, 11, 12];

fprintf('原始矩阵 A =\n');
disp(A);

% Givens QR分解
[Q, R] = givens_qr_decomp(A);

fprintf('Q矩阵 =\n');
disp(Q);
fprintf('R矩阵 =\n');
disp(R);

% 验证分解
QR_product = Q * R;
fprintf('验证 Q*R - A 的误差: %e\n', norm(QR_product - A, 'fro'));

% 验证Q的正交性
QtQ = Q' * Q;
fprintf('Q的正交性 ||Q^T*Q - I||: %e\n', norm(QtQ - eye(size(QtQ)), 'fro'));

%% 4. 稀疏矩阵的Givens QR分解
fprintf('\n4. 稀疏矩阵的Givens QR分解\n');
fprintf('Givens旋转特别适合稀疏矩阵，因为它只影响两行\n\n');

% 创建稀疏三对角矩阵
n = 6;
A_sparse = diag(2*ones(n, 1)) + diag(-ones(n-1, 1), 1) + diag(-ones(n-1, 1), -1);
fprintf('稀疏三对角矩阵 A =\n');
disp(full(A_sparse));
fprintf('非零元素数: %d / %d (%.1f%%)\n', nnz(A_sparse), numel(A_sparse), ...
        100 * nnz(A_sparse) / numel(A_sparse));

% Givens QR分解
[Q_sparse, R_sparse] = givens_qr_decomp(A_sparse);

fprintf('R矩阵的非零元素数: %d / %d (%.1f%%)\n', nnz(R_sparse), numel(R_sparse), ...
        100 * nnz(R_sparse) / numel(R_sparse));

% 验证
fprintf('稀疏矩阵QR分解误差: %e\n', norm(Q_sparse * R_sparse - A_sparse, 'fro'));

%% 5. 与Householder QR的比较
fprintf('\n5. 与Householder QR的比较\n');
fprintf('比较Givens和Householder方法的特点\n\n');

% 使用相同的测试矩阵
A_test = randn(8, 5);

% Givens QR
tic;
[Q_givens, R_givens] = givens_qr_decomp(A_test);
time_givens = toc;

% Householder QR (简化版本)
tic;
[Q_house, R_house] = householder_qr_simple(A_test);
time_house = toc;

fprintf('计算时间比较:\n');
fprintf('Givens QR: %.6f 秒\n', time_givens);
fprintf('Householder QR: %.6f 秒\n', time_house);
fprintf('时间比: %.2f\n', time_givens / time_house);

% 精度比较
error_givens = norm(Q_givens * R_givens - A_test, 'fro');
error_house = norm(Q_house * R_house - A_test, 'fro');

fprintf('分解精度比较:\n');
fprintf('Givens误差: %e\n', error_givens);
fprintf('Householder误差: %e\n', error_house);

%% 6. 快速Givens旋转
fprintf('\n6. 快速Givens旋转\n');
fprintf('避免计算三角函数的快速Givens变换\n\n');

% 演示快速Givens旋转
a = 3; b = 4;
fprintf('原始元素: a = %g, b = %g\n', a, b);

% 标准Givens旋转
[c_std, s_std] = givens_rotation(a, b);
result_std = [c_std, s_std; -s_std, c_std] * [a; b];

% 快速Givens旋转
[gamma, sigma, result_fast] = fast_givens_rotation(a, b);
fprintf('快速Givens参数: γ = %.6f, σ = %.6f\n', gamma, sigma);
fprintf('标准结果: [%.6f; %.6f]\n', result_std);
fprintf('快速结果: [%.6f; %.6f]\n', result_fast);
fprintf('结果差异: %e\n', norm(result_std - result_fast));

%% 7. 应用：上Hessenberg矩阵的QR分解
fprintf('\n7. 应用：上Hessenberg矩阵的QR分解\n');
fprintf('上Hessenberg矩阵的QR分解只需要n-1次Givens旋转\n\n');

% 创建上Hessenberg矩阵
n = 5;
H = triu(randn(n, n), -1);  % 上Hessenberg形式
fprintf('上Hessenberg矩阵 H =\n');
disp(H);

% 专门的Hessenberg QR分解
[Q_hess, R_hess] = hessenberg_qr(H);

fprintf('Hessenberg QR分解误差: %e\n', norm(Q_hess * R_hess - H, 'fro'));
fprintf('所需的Givens旋转次数: %d\n', n-1);

%% 8. 并行化潜力
fprintf('\n8. 并行化潜力\n');
fprintf('Givens旋转的并行化特性\n\n');

% 演示可以并行执行的Givens旋转
A_para = randn(6, 4);
fprintf('演示矩阵 A (%d×%d):\n', size(A_para));
disp(A_para);

% 分析并行化模式
[parallel_stages, total_rotations] = analyze_givens_parallelization(A_para);
fprintf('总Givens旋转数: %d\n', total_rotations);
fprintf('并行阶段数: %d\n', length(parallel_stages));

for stage = 1:length(parallel_stages)
    fprintf('阶段 %d: %d 个并行旋转\n', stage, parallel_stages(stage));
end

%% 9. 性能测试
fprintf('\n9. 性能测试\n');
fprintf('不同规模矩阵的Givens QR性能\n\n');

sizes = [50, 100, 200];
fprintf('矩阵规模    Givens QR    运算量估计\n');
fprintf('--------    ---------    ----------\n');

for i = 1:length(sizes)
    m = sizes(i);
    n = min(m, 50);
    A_perf = randn(m, n);
    
    tic;
    [Q_perf, R_perf] = givens_qr_decomp(A_perf);
    time_perf = toc;
    
    % 理论运算量：约 2mn² - (2/3)n³ 次乘法
    flops_estimate = 2*m*n^2 - (2/3)*n^3;
    
    fprintf('%6d      %7.4f s    %8.0f\n', m, time_perf, flops_estimate);
end

end

function [c, s] = givens_rotation(a, b)
% GIVENS_ROTATION 计算Givens旋转参数
% 
% 输入: a, b - 要消除的元素
% 输出: c, s - 旋转参数，使得 [c s; -s c] * [a; b] = [r; 0]

if b == 0
    c = 1;
    s = 0;
elseif abs(b) > abs(a)
    tau = -a / b;
    s = 1 / sqrt(1 + tau^2);
    c = s * tau;
else
    tau = -b / a;
    c = 1 / sqrt(1 + tau^2);
    s = c * tau;
end

end

function [Q, R] = givens_qr_decomp(A)
% GIVENS_QR_DECOMP Givens旋转QR分解
% 
% 输入: A - m×n 矩阵
% 输出: Q - m×m 正交矩阵, R - m×n 上三角矩阵

[m, n] = size(A);
Q = eye(m);
R = A;

% 逐列消除下三角元素
for j = 1:n
    for i = m:-1:j+1
        if abs(R(i, j)) > eps
            % 计算Givens旋转参数
            [c, s] = givens_rotation(R(i-1, j), R(i, j));
            
            % 构造Givens矩阵
            G = eye(m);
            G([i-1, i], [i-1, i]) = [c, s; -s, c];
            
            % 应用旋转
            R = G * R;
            Q = Q * G';
        end
    end
end

end

function [gamma, sigma, result] = fast_givens_rotation(a, b)
% FAST_GIVENS_ROTATION 快速Givens旋转
% 
% 避免计算三角函数的快速算法

if b == 0
    gamma = 1;
    sigma = 0;
    result = [a; 0];
elseif abs(b) > abs(a)
    tau = a / b;
    sigma = 1 / sqrt(1 + tau^2);
    gamma = sigma * tau;
    result = [b / sigma; 0];
else
    tau = b / a;
    gamma = 1 / sqrt(1 + tau^2);
    sigma = gamma * tau;
    result = [a / gamma; 0];
end

end

function [Q, R] = hessenberg_qr(H)
% HESSENBERG_QR 上Hessenberg矩阵的QR分解
% 
% 专门针对上Hessenberg矩阵的高效QR分解

n = size(H, 1);
Q = eye(n);
R = H;

% 只需要消除下对角线元素
for i = 2:n
    if abs(R(i, i-1)) > eps
        [c, s] = givens_rotation(R(i-1, i-1), R(i, i-1));
        
        G = eye(n);
        G([i-1, i], [i-1, i]) = [c, s; -s, c];
        
        R = G * R;
        Q = Q * G';
    end
end

end

function [stages, total] = analyze_givens_parallelization(A)
% ANALYZE_GIVENS_PARALLELIZATION 分析Givens旋转的并行化
% 
% 分析QR分解中可以并行执行的Givens旋转

[m, n] = size(A);
stages = [];
total = 0;

for j = 1:n
    rotations_in_column = m - j;
    if rotations_in_column > 0
        % 在每列中，旋转可以从底部开始并行执行
        stage_count = 0;
        remaining = rotations_in_column;
        
        while remaining > 0
            parallel_rotations = floor(remaining / 2);
            if parallel_rotations == 0
                parallel_rotations = 1;
            end
            stages(end + 1) = parallel_rotations;
            stage_count = stage_count + 1;
            remaining = remaining - parallel_rotations;
        end
        
        total = total + rotations_in_column;
    end
end

end

function [Q, R] = householder_qr_simple(A)
% HOUSEHOLDER_QR_SIMPLE 简化的Householder QR分解
% 
% 用于与Givens方法比较的简化版本

[m, n] = size(A);
Q = eye(m);
R = A;

for k = 1:min(m-1, n)
    x = R(k:m, k);
    
    % 计算Householder向量
    v = x;
    v(1) = v(1) + sign(x(1)) * norm(x);
    v = v / norm(v);
    
    % 应用Householder变换
    H = eye(m-k+1) - 2 * (v * v');
    R(k:m, k:n) = H * R(k:m, k:n);
    Q(k:m, :) = H * Q(k:m, :);
end

Q = Q';

end